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I.  INTRODUCTION 


High-current  pulsed  electron  beam  accelerator  concepts  are  being 

investigated  in  our  laboratory.  One  particularly  interesting  concept1 
is  based  on  accelerating  sections  containing  charged  internally-switched 
elements.  By  closing  the  various  switches  in  the  proper  synchronization, 
electromagnetic  pulses  are  produced  that  would  accelerate  pulsed  beams . 

2 3 

This  concept  has  recently  been  analyzed  ’ by  considering  the  cavities 
in  which  the  electromagnetic  pulses  travel  as  transmission  lines  and  ap- 
plying principal  mode  theory.  It  was  soon  realized  that  geometries  that 
correspond  to  constant- impedance  transmission  lines  are  advantageous  for 
producing  uniformly  accelerated  beams.  Pulse  distortions  at  the  line 
discontinuities  are  being  investigated  by  employing  equivalent  circuit 
2 4-7 

representation  ’ of  the  discontinuities.  Although  such  an  analysis 
provides  a closed  form  representation  to  the  time  dependence  of  the  re- 
sulting accelerating  pulse,  it  leaves  much  to  be  desired.  Firstly,  as 
the  frequency  is  increased  the  lumped  parameters  become  frequency  depen- 
dent and  above  a critical  frequency  the  entire  equivalent  circuit  repre- 
sentation breaks  down.  Secondly,  complications  arise  when  discontinuities 
are  close  to  each  other.  And  thirdly  the  value  of  the  equivalent  para- 
meters is  not  available  for  many  interesting  geometries.  Another  feature 
not  addressed  by  the  transmission  line  analysis  is  the  coupling  of  the 
beam  to  the  accelerating  structure. 

In  order  to  complement  the  transmission  line  analysis  and  address 
the  above  issues,  computer  simulations  with  electron  beams  have  been  per- 
formed. A two-dimensional  axi-symmetric  code  was  obtained  from  the  Naval 
Research  Laboratory  that  directly  integrates  Maxwell's  equations  and  then 
accelerates  charged  particles  with  the  resulting  fields.  Additions  and 
modifications  were  made  to  the  code  so  that  it  would  address  some  of  the 
questions  raised  by  the  transmission- line  analysis. 
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Section  II  will  present  the  formulation  and  limitations  of  the  code, 
and  Section  III  will  describe  how  the  code  is  applied  to  the  current 
problem.  The  remainder  of  this  report  is  devoted  to  the  results,  dis- 
cussion and  direction  for  future  work.  A program  users  guide  is  provided 
in  the  Appendix. 


II.  FORMULATION 


A.  Electromagnetic  Fields 

The  propagation  of  electromagnetic  fields  is  governed  by  Maxwell's 

8 

equations.  In  their  usual  form  the  equations  are  written  as  (we  use 
Gaussian  units  since  these  are  the  units  for  which  the  code  was  written) 


V*D  = 4irp 


V'3  = 0 


» x fl  . 4L } . I |I 

c c at 


and 


V x 


c at 


= o 


where 


(1) 

(2) 

(3) 

(4) 


D is  the  electric  displacement, 

^ is  the  electric  field, 

£ is  the  magnetic  flux  density, 
ft  is  the  magnetic  field, 

J is  the  current  density, 

p is  the  charge  density, 

and  c is  the  velocity  of  light. 

~0  " ' 

Classical  Electrodynamics , J.D.  Jackson,  p.  178,  J,  W%le y,  1962. 
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In  vacuo 


5 = t 

and  (5) 

t = H . 


We  will  consider  only  the  in  vacuo  case  for  this  report. 

If  the  initial  values  of  the  fields  are  given  and  knowledge  of 
their  time  evolution  is  desired,  it  is  convenient  to  write  eqs.  (3  and 
4)  using  eqs.  (5)  as 

|f  = c V x £ - 471?  (6) 


and 


3$ 

3t 


c V x E. 


(7) 


Application  of  the  continuity  equation 


7.J  + |£.  = 0 (8) 


to  eqs.  (.6  and  7)  will  show  that  if  eqs.  (1  and  2)  are  initially  obeyed 
they  will  continue  to  be  obeyed  when  subjected  to  the  time  evolution 
prescribed  by  eqs.  (6  and  7);  thus  a time  integration  algorithm  need  only 
consider  eqs.  (6  and  7). 

A staggered- leapfrog  scheme  similar  to  that  of  Boris9  is  used  to 
advance  the  E and  B fields.  The  space  grids  used  are  shown  in  Fig.  1. 

For  our  present  purposes  only  the  transverse  magnetic  (TM)  modes  are 
considered.  This  restriction  combined  with  the  two-dimensional  axi- 
symmetric  nature  of  the  code  implies  that  only  the  E^,  ER,  B0 , JR,  and 

J field  components  enter  into  the  calculations.  The  J field  component 
Z K 

J.P.  Boris,  "Relativistic  Plasma  Simulation-Optimization  of  a Hybrid 
Code , " in  the  Proceedings  of  the  Fourth  Conference  on  Numerical  Simula- 
tion of  Plasmas , Ed.  by  J.P.  Boris  and  Rama  C.  Shanny , pp.  3-67,  Naval 
Research  Laboratory,  1970. 
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Figure  1.  Space  grid  for  the  finite  difference  Maxwell  equations 
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is  suppressed  because  allowing  the  particles  a radial  degree  of  freedom 
would  lead  to  beam  blow-up  if  no  stabilizing  field  is  present.  Such 

a field  of  course  introduces  non-'IN  modes.  Later  we  will  justify 

and  explain  the  motivation  for  limiting  the  calculations  to  the  TM  modes. 
The  fields  are  evolved  with  the  following  equations: 


Ez  (Z  - AZ/2,  R,  T + AT/2)  = E^Z  - AZ/2,  R,  T - AT/2)  + 


CAT 

AR 


2R  + AR 
2R 


Bq  (Z  - AZ/2,  R + AR/2,  T) 
0 


2R  - AR 
2R 


B^  (Z  - AZ/2,  R 


AR/2 


, T) 


(9) 


- AT  Jz  (Z  - AZ/2,  R,  T)  , 

E (Z,  R - AR/2,  T + AT/2)  = ED  (Z,  R - AR/2,  T - AT/2) 

K K 

- ~ - [B0  (Z  + AZ/2,  R - AR/2,  T)  (10) 

- B0  (Z  - AZ/2,  R - AR/2,  T) ] , 


and 


Bq  (Z  - AZ/2,  R - AR/2,  T + AT)  = B (Z  - AZ/2,  R - AR/2,  T) 

0 u 

+ ^ [Ez  (Z  - AZ/2,  R,  T + AT/2)  - Ez  (Z  - AZ/2,  R - AR,  T ♦ 
- ^ [Er  (Z,  R - AR/2,  T + AT/2)  - ER  (Z  - AZ,  R - AR/2,  T + 
On  the  axis  (R=0)  eq.  (9)  is  replaced  by 


AT/2)]  (11) 
AT/2)]  . 
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E (Z  - AZ/2,  0,  T + A/2)  = E (Z  - AZ/2,  0,  T - AT/2)  + 

L * Lt 

Bq  (Z  - AZ/2,  AR/2,  T)  (12) 

- 4ttAT  Jz  (Z  - Z/2,  0,  T)  . 


This  equation  follows  from  an  application  of  Stokes'  theorem.  Also  on 

the  axis  B and  ED  are  set  to  zero.  This  is  done  in  practice  by  setting 
0 K 

Er  (Z,  - AR/2,  T + AT/2)  = - ER(Z,  AR/2,  T + AT/2)  (13) 


and 


Bq  (Z  - AZ/2,  - AR/2,  T)  = - BQ  (Z  - AZ/2,  AR/2,  T)  . (14) 

0 t) 

The  code  at  present  allows  only  those  boundaries  which  lie  along 
the  Z or  R directions.  Perfectly  conducting  boundaries  are  assumed; 
thus  the  tangential  components  of  E and  the  normal  component  of  B are 
zero  on  the  boundaries.  With  our  restricted  number  of  field  components, 
only  the  boundary  conditions  on  E enter  into  consideration.  The  code 
sets  the  E field  to  zero  within  the  conductors  by  means  of  rectangular 
(in  an  R,  Z plane)  masks.  Two  masks  are  required,  one  for  E and  one 

for  E^,  because  of  the  different  boundary  conditions  and  grid  point  lo- 
cations of  these  field  components. 

The  solution  of  the  finite  difference  equations  (9-14)  are  not  the 
exact  solutions  to  Maxwell's  equations.  An  effective  way  to  judge  how 
good  are  the  finite  difference  solutions  is  to  examine  their  dispersion 
relation.  First  we  will  examine  the  dispersion  relations  for  the  dif- 
ferential equations  so  that  we  will  have  a basis  for  comparison.  From 
Maxwell's  equation  it  is  a simple  matter  to  show  that,  in  a charge  free 
region,  the  electric  field  obeys  the  wave  equation, 


(15) 


For  our  geometry  we  obtain  the  following  eouations  for  the  field  com- 

* 10 
ponents; 

Methods  of  Theoretical  Physios,  P.M.  Morse  and  H.  Feshbach,  p.  116, 
McGraw-Hill,  1953. 
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i 


ii(R\  s + 

R 3R  1 3R  R2  3Z2 


c2  3t‘ 


I _2  fR  Ziw 

R 3R  L 3RJ 


i »\ 
c2  3,2 


Now  assume  the  field  components  have  the  form  Rn  exp  [i(cot  + k R + 
Z)]  . Equations  (16)  and  (17)  then  become  respectively:  R 


2 ->  ikn  . 2 

\ = c2  [1  - (2n  + 1) 
k k R k2R2 


^ = c2  [1  - (2n  + 1)  - ~~T~2  J * (19) 

k kV  kV 

where  k2  = k2  + k2  . 

If  kR  = 0,  then  any  propagating  part  of  the  field  must  have  E z = 0; 

and  it  is  convenient  to  choose  n = - 1 so  that  eq.  (18)  becomes  nondis- 
2 2 2 

persive,  a>  /k  = c . (The  choice  n = 1 would  violate  the  equation 
V-£  = 0 that  results  from  our  charge-free  assumption.)  If  k ^ 0,  then 
for  stability  reasons  to  keep  the  dispersion  relations  real,  n must  be 
- j.  In  this  case  eqs.  (18)  and  (19)  become 

2 ? , 

^2  = c (1  + — TT  ) (2°) 

k 4k  R 


•1  „ c 2 CI  - _i— 
2 *•  2 2 
k 4k  R 


Now  back  to  the  difference  equations.  Just  as  with  the  differential 
equations,  it  is  most  convenient  to  derive  the  dispersion  relations  from 
a wave  equation.  We  now  derive  the  finite  difference  form  of  the  wave 
equation  directly  from  the  finite  difference  form  of  Maxwell's  equations. 
The  notation  may  be  simplified  by  introducing  the  following  operators: 
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3t  X (Z,  R,  T)  = [X(Z,  R,  T + AT/2)  - X(Z,  R,  T - AT/2)]/AT,  (22) 

3Z  X (R,  Z,  T)  = [X(Z  + AZ/2,  R,  T)  - X(Z  - AZ/2,  R,  T)]/AZ,  (23) 

3r  X(R,  Z,  T)  = [X(Z,  R 1-  AR/2,T)  - X(Z,  R - AR/2,T)]/AR,  (24) 

and 

ad 

3'r  X(R,Z,T)  = [(1  ♦ ^)X(Z,R  + AR/2.T) 

- (1  - §f)X(Z,R  - AR/2,T)]/AR,  (25) 

where  X is  any  field  component. 

We  note  that  all  pairs  of  these  operators  commute  except  3 and  3R. 
Equations  (9-11)  can  now  be  written  compactly  as 

aT  Ez  ■ c SR  Be  - 4*JZ' 

3r  er  ■ - 43z  B0  • C26) 

and 


3t  Be  ■ cI3r  Ez  - az  V • 

In  a charge-free  region  V*£  = 0;  the  finite  difference  form  of  this  equa 
tion  is 

3RER  ♦ 3z  Ez  ■ 0 • (27) 

From  eqs.  (26,27)  follow  in  a charge-free  region 


~~2  3R  SR  ER  * 3Z  ER 

C 


(28) 
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and 


2 


SR  SE  EZ  * 3Z  EZ 


These  last  two  equations  expand  to  read 


(29) 


~2  [Er(Z,R,T  + AT)  - 2 Er(Z,R,T)  + ER  (Z,R,  T - AT)]/(AT)2  = 

1(1  * #Tffi)ERcz.R  * T)  * - zir^AR  - 2>er<z^'t> 

AD  9 

+ (1  ' 2R~-  AR)ER(Z>R  ’ AR’  T)J/(AR)  + 


[Er(Z  ♦ AZ.R.T)  - 2 Er(Z,R,T)  + ER(Z  - AZ,R,T)]/(AZ) 2 (30) 


and 

~2  tEz(z»  R.T  + AT)  - 2EZ(Z,R,T)  + EZ(Z,R,T  - AT)]/(AT)2  = 

AR 

1(1  + 2R°  EZ(Z*R  + AR*T)  _ 2 EZ(Z’R>T)  (31) 

♦ (!  - §§0  EZ(Z,R  - AR,T) ]/ (AR) 2 

+(EZ(Z  + AZ,R,T)  - 2 EZ(Z,R,T)  + Ez  (Z  - AZ.R.T) ]/ (AZ) 2 . 

Again  we  assume  that  the  field  components  have  the  form  Rn  exp 
[i(u>t  + kR  R + kzZ)].  As  before  when  kR  = 0,  we  take  n = -1;  and  the 

resulting  form  for  ER  inserted  into  eq.  (30)  yields 
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When  k / 0,  we  take  as  before  n = - 
K 

into  eq.  (30)  we  obtain 


1 

V 


and  when  this  form  is  inserted 


1 / iwAT  . -iwAA  _ 1 li 

(CAT)2  \ 6 / (AZ) 2 \ 

1 T/,  AR  \ /.  AR\ 

Vr)2  [(  2R*al!)(  ' 


ik_  AZ  + -ik?  AZ 
Z - l + e Z 


) 


1/2  eikR  AR  - 2 ♦ -** 


2R  + AR 


(33) 


AR 


2R  - AR 


0-1) 1/2  *-ikR“  • 


The  coefficient  of  the  exponentials  are  then  expanded  in  power  series  of 
AR/R  and  the  resulting  series  are  multiplied  term  by  term.  Proceeding  in 
a straightforward  manner  after  collecting  the  common  powers  of  AR/R  yields 
the  following  dispersion  relation  for  ER 


where  is  the  binomial  coefficient.  Proceeding  in  the  same  manner 
from  eq.  (31)  we  obtain  the  dispersion  relation  for  E^  as  follows. 


+ 


cos(k„R) 

i sin(kR 
2 

16R 


(35) 


(34) 
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In  the  limit  AT,  AR,  and  AZ  -*■  0,  eqs.  (34-35)  reduce  to  eqs.  (20-21) 
as  expected.  Now  if  we  do  not  want  the  numeric  solution  to  have  non- 
physical exponential  time  growth,  u must  be  real.  The  last  term  in  eqs. 
(34-35)  will  contribute  an  imaginary  part  to  to  if  kR  and  k z are  real, 

but  this  contribution  is  only  significant  near  the  axis.  As  the  pulse 
leaves  ♦'his  region  (remember  now  kR  / 0)  , AR/R  -*■  0 and  this  last  term 

vanishes.  What  these  terms  mean  is  that  we  have  not  chosen  exactly  cor- 
rect radial  dependences  near  the  axis  for  the  finite  difference  equations. 

-1/2 

The  correct  radial  dependence  is  not  R near  the  axis  but  does  approach 


A most  serious  contribution  to  the  imaginary  part  of  io  would  come 

2 

from  having  the  right  side  of  eqs.  (34-35)  larger  than  l/(cAT)  since 

then  sin  will  become  larger  than  unity.  Thus  to  insure  stability 

for  all  values  of  k and  k the  time  steps  must  be  limited  by  the  relation 
K Z 


— + Q , (36) 

(cAT)  (AR)  (AT) 

where  Q consists  of  the  terms  containing  the  summations  in  eqs.  (34-35). 
Usually  we  may  neglect  Q since  it  is  generally  negative. 

An  important  consequence  of  the  difference  between  the  exact  dis- 
persion relations  and  the  finite  difference  dispersion  relations  is  that 
although  the  lower  frequency  finite  difference  wave  propagates  at  essen- 
tially the  velocity  of  light,  the  higher  frequency  finite  difference 

9 

waves  propagate  more  slowly.  Thus  square  pulses  would  have  their  rise- 
times  degraded  and  would  develop  overshoots  and  ringing.  These  annoying 
numeric  artifacts  may  be  minimized  by  smoothing  out  the  rise  and  fall  of 
an  injected  square  pulse  and  taking  a fine  mesh  if  the  required  computer 
resources  are  available.  It  has  been  determined  empirically  that  the 
numeric  overshoots  reach  a maximum  amplitude  above  the  square  pulse  of 
(4/ir-l)  times  the  amplitude  of  the  square  pulse  if  no  reflections  from 
boundaries  occur.  Reflections  increase  the  amount  of  overshoot.  We  note 
in  passing  that  4/ir  is  the  amplitude  of  the  lowest  frequency  component 
of  a unit  amplitude  square  pulse. 

We  now  consider  the  formulation  for  the  solution  of  Poisson's  equa- 
tion 


(37) 
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where  4>  is  the  electric  potential.  This  equation  arises  when  we  would 
like  to  know  the  potential  and  electric  field  components  due  to  the  ini- 
tial charging  up  of  electrodes  within  the  cavities.  Usually  p = 0;  how- 
ever, the  general  case  can  be  handled  with  a minimum  amount  of  additional 
effort.  Reference  11  contains  a good  discussion  of  relaxation  and  over- 
relaxation techniques  for  solving  the  finite  difference  equation  corres- 
ponding to  eq.  (37). 

The  potential  is  specified  on  all  points  of  a closed  boundary  and 
the  potential  at  the  interior  points  for  our  geometry  are  obtained  from 
successive  application  of  the  following  formulae 

$(n+1)(Z,R)  = *(n)(Z,R)  + 

, n ( 1 fVn)(Z  + AZ,R)  + <fr(n)(Z  - AZ,R) 

(2/(AZr  + 2/(ARr  L (AZ) 

+ R * AR/2  tp  (Z,  R + AR)  + R-  AR/-2 
R(AR)  R (AR) ^ 

- /n)(Z,R)|. 

/ 

On  the  axis  we  use 


4>(n)(Z,R  - AR)  - 4irp  (Z ,R)1  (38) 


/n+1)(z,o)  = <(> (n^  (Z, 0)  + 


5f 


2/ (A  Z) z + 4/ (AR) ' 


1 


<t(n)(z 


+ AZ , 0)  + $(n)(Z  - AZ , 0) 


(AZ) 


+ -J—  $(n)(Z,AR)  - 4tt p (Z , 0)  I - 4>(n)(Z,0)l  . (39) 

(AR)2  J ) 

An  overrelaxation  factor  of  6 = 1.7  was  found  to  be  satisfactory  for  our 
application.  We  iterated  with  eqs.  (38-39)  until  <J>  changed  by  less  than 

_7 

a specified  factor  of  usually  10  times  the  maximum  <p.  Having  obtained 
<j> , the  electric  field  components  are  computed  from 


Space-Charge  Flou),  Peter  T.  Kirstein , George  S.  Kino , and  William  E. 
Walters,  Chap.  9,  McGraw-Hill,  1967.  Some  typographical  errors  exist, 
in  this  chapter.  Of  interest  to  us  are  the  following:  the  factor  h 
should  be  removed  from  the  denominator  before  the  brackets  in  e<^.(3.1) 
and  the  last  terms  of  eqs.  (3.2,  3.4)  should  be  multiplied  by  n. 
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(40) 


Ed(Z,R  - AR/2)  = ^ ^R— 


and 


EZ(Z  - AZ/2.R)  = ^ 


(41) 


B.  Particle  Pushing 

The  particle  motion  is  prescribed  by  relativistic  equations  of  mo- 
tion. All  of  the  particle  motion  is  constrained  to  be  in  the  Z direction 
for  reasons  to  be  given  later.  Under  this  constraint  the  E_  field  com- 

Lt 

ponent  is  the  only  computed  field  component  that  provides  an  acceleration 
in  the  Z direction.  The  value  of  E^  at  a particle  position  is  calculated 

by  interpolating  the  values  of  E^  from  the  4 nearest  mesh  points.  The 

particle  momentum  per  unit  mass  is  advanced  by 

P (T > I ) = P (T  - AT, I)  ♦ ^ Eint  (T  - AT/2)  AT  (42) 

where  I is  the  particle  index, 

Q is  the  particle  charge, 

M is  the  particle  mass,  and 
£1^  is  the  interpolated  E^ . 

The  particle  positions  are  then  advanced  by 

Z (T  + AT/2,1)  = Z (T  - AT/2,1)  + P(T,I.2_AT  . (43) 

l . [EHjii]2 

Each  particle  contributes  to  the  and  p fields  by  an  inverse  interpola- 
tion procedure;  that  is,  the  particle's  charge  and  current  are  apportioned 
to  the  nearest  corresponding  4 grid  points  with  the  amount  of  apportion- 
ment depending  on  how  close  the  particle  is  to  a grid  line. 

A soft  injection  of  the  particles  is  employed.  This  means  that  the 
particles  are  "created"  with  their  initial  energy  one-half  of  grid  spac- 
ing before  the  beginning  of  the  field  grids.  Until  these  so-called  in- 
active particles  reach  the  beginning  of  the  field  lattice  these  particles 
do  not  respond  to  the  fields- -that  is,  they  are  not  accelerated;  however 
they  still  contribute  to  the  p and  fields.  Likewise  upon  leaving  the 

field  lattice,  the  particle  still  contribute  to  the  p and  J fields  but 
are  not  accelerated  for  one  half  lattice  spacing  at  the  end  of  the  field 
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lattice.  We  will  say  more  about  these  transition  regions  in  the  dis- 
cussion section. 


III.  IMPLEMENTATION 

A cross-section  al  view  of  the  cylindrical ly- symmetric  cavity  to  be 
simulated  is  shown  in  Figure  2.  Initially  the  cylinder  within  the  cav- 
ity is  placed  at  a specified  voltage  and  all  of  the  other  surfaces  are 
held  at  ground  potential.  The  switch  is  considered  nonexistent  at  this 
point  in  the  calculation.  Then  the  Poisson's  equation  solver  subroutine 
calculates  the  resulting  potential  and  electric  field  components.  At 
this  point  we  of  course  have  VX  £ = 0.  Then  closing  of  the  switches  is 
simulated  by  setting  the  electric  field  components  to  zero  in  the  region 
of  the  switches.  This  causes  non-zero  values  of  V x £ to  appear  which 
act  as  driving  terms  for  eqs.  (7)  and  (11)  and  initiate  the  propagating 
electromagnetic  pulse.  As  mentioned  before  and  as  experience  has  shown, 
immediately  setting  the  electric  field  to  zero  gives  rise  to  undesirable 
overshoots  as  the  pulse  propagates.  The  overshoot  may  be  reduced  al- 
though not  completely  eliminated  by  setting  the  field  to  zero  slowly. 

In  practice  we  do  this  by  setting 


£(T)  = exp(-T/T0), 


(44) 


where  EQ  is  the  initial  electric  field  and  TQ  is  a decay  time  constant. 

Figure  3 shows  the  rise  of  the  overshoot  as  a function  of  time  step  for 
Tq  = 15  AT  as  calculated  from  a one-dimensional  analog  of  eqs.  (9  and 

11).  In  the  way  of  contrast,  for  TQ  ■*  0 the  overshoot  approaches  25% 
after  only  120  time  steps. 

Initial  calculations  are  performed  without  injection  of  a particle 
beam.  The  progress  of  the  electromagnetic  pulse  propagation  is  diag- 
nosed by  periodic  listings  and  plotting  of  all  or  selected  field  compon- 
ents. In  addition,  the  open-circuit  voltage  across  the  accelerating  gap 
is  monitored  as  well  as  the  electric,  magnetic  and  total  energies.  The 
open  circuit  voltage  as  a function  of  time  is  stored  on  disk  for  future 
processing  and  for  plotting. 

After  the  open-circuit  calculations  have  been  performed  and  the  pre- 
cise timing  of  the  arrival  of  the  accelerating  pulse  is  known,  calcula- 
tions are  performed  with  beam  injection.  The  timing  of  the  beam  injec- 
tion with  respect  to  the  accelerating  pulse  is  critical  if  overshoots 
and  ringing  in  the  accelerated  beam  energy  are  unwanted.  The  beam  cur- 
rent is  tailored  to  have  a slow  exponential  rise  and  fall  so  that  large 
dl/dT  voltages  are  not  generated  as  the  beam  passes  the  gap.  After  the 
beam  passes  the  accelerating  gap,  its  energy  is  calculated  at  selected 
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ers 


time  intervals  and  radial  positions,  and  this  information  is  stored  on 
disk  for  subsequent  handling.  Both  solid  and  hollow  beams  can  be  em- 
ployed. 


IV.  RESULTS 

Figure  2 shows  the  geometry  used  for  most  of  the  calculations.  Some 
small  modifications  were  however  made  to  test  their  effect  and  improve 
the  cavity  performance  as  will  be  mentioned  later.  The  length  of  the 
cavity  was  chosen  to  produce  an  approximately  10  nsec  accelerating  pulse, 
and  diameters  were  chosen  to  provide  reasonable  proportions  and  to  pro- 
vide proper  impedance  ratios.  The  outer  and  inner  portions  of  the  cav- 
ity have  characteristic  impedances  of  20.2  and  61.8  ohms.  These  are  in 
a ratio  of  3.06  to  1 as  compared  to  a 3 to  1 ratio  that  the  transmis- 
(2  3) 

sion  ’ line  analysis  concludes  is  one  desirable  possibility.  Figure 
4 shows  the  open-circuit  voltage  expected  from  the  principal  mode  trans- 
mission line  analysis  when  the  effects  of  discontinuity  at  the  right  end 
of  the  cavity  and  at  the  gap  are  neglected.  With  our  simulation  code 
the  open  circuit  voltage  shown  in  Figure  5 is  obtained  when  the  initial 
charge  on  the  middle  cylinder  is  1000  statvolts  O 300  kV) . Other  parame- 
ters of  this  simulation  are  AT  = 10'11  sec,  AR  = AZ  = 1 cm,  and  T = 

-10  0 
1.5  x 10  sec.  Some  of  the  lower  amplitude  oscillation  is  due  to  the 

previously  mentioned  numeric  dispersion,  but  the  larger  amplitude  oscilla- 
tions at  0,  10,  20-22,  30-32,  and  beyond  40  ns  are  physical  effects  due 
to  impedance  discontinuities  and  excitation  of  higher  order  modes  from 
the  switch  closing.  Figure  6 shows  an  intensity  field  of  B 4 ns  after 

the  switch  closing.  A bipolar  pulse  due  to  the  switch  closing  is  clearly 
seen  in  the  inner  transmission  line.  Such  a pulse  is,  of  course,  not  in- 
cluded in  the  idealized  transmission  line  analysis  and  gives  rise  to  some 
of  the  structure  in  the  open  circuit  voltage. 

Of  particular  interest  is  how  much  the  risetime  of  the  output  volt- 
age pulse  is  degraded.  The  most  useful  single  pulse  for  acceleration  is 
that  between  about  10  and  20  ns.  The  polarity  change  at  about  10  ns  pro- 
ceeds from  10%  to  90%  of  the  excursion  in  0.6  ns.  This  corresponds  to 
an  exponential  rise  time  of  0.27  ns.  The  excursion  around  20  ns  takes 
somewhat  longer.  To  compare  with  equivalent  circuit  models  of  discon- 
tinuities it  is  also  of  interest  to  study  the  risetime  degradation  of 
the  pulse  passing  just  one  discontinuity,  for  example,  the  right-hand 
end  of  the  cavity.  Such  effects  cannot  be  obtained  with  the  geometry  of 
Fig.  2 since  there  are  two  pulses  propagating  in  opposing  directions  that 
pass  the  discontinuity  at  the  same  time.  Figure  7 illustrates  the  modi- 
fication of  the  geometry  employed  to  overcome  this  difficulty.  The 
baffles  in  the  inner  transmission  line  provide  a tortuous  path  that  slow 
the  pulse  traveling  through  it  so  that  the  pulse  in  the  outer  transmis- 
sion line  may  pass  around  the  discontinuity  before  the  other  pulse  ar- 
rives in  the  vicinity.  An  intensity  plot  of  B^  with  this  geometry  at  a 
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Figure  6.  B 4 ns  after  switch  closing.  Every  other  point  is  plotted  in  the  Z direction. 
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time  of  (■>  ns  is  shown  in  figure  8.  The  risetime  of  the  simulated  voltage 
pulse  just  before  going  around  the  discontinuity  is  0.16  ns  and  just  after 
it  is  0.27  ns.  Unfolding  of  the  "risetime"  of  the  discontinuity  from 
these  numbers  is  somewhat  ambigious  since  a circuit  containing  several 
subcircuits  each  having  its  own  risetime  has  an  effective  rise  that  is 
not  a simple  function  of  the  subcircuit  risetimes.  The  total  circuit 
risetime  is  less  than  the  sum  of  the  individual  risetimes  and  approaches 
the  square  root  of  the  sum  of  the  squares  of  the  individual  risetimes 

1 2 

when  there  are  many  subcircuits.  A linear  subtraction  of  the  above 
risetimes  gives  0.11  ns  while  subtraction  in  quadrature  yields  0.22  ns 
for  the  risetime  of  the  discontinuity.  The  actual  value  is  probably  much 
closer  to  the  second  figure  but  somewhat  less  than  it. 

After  having  simulated  the  electromagnetic  pulse  in  the  cavity,  we 
include  the  injection  of  an  electron  beam  from  the  left-hand  side  of  the 
inner  cylinder.  Solid  and  hollow  beams  were  tried,  the  latter  givinc 

2 

slightly  better  results.  The  beam  energy  gained  in  units  of  mgc  (511 
keV)  is  shown  as  a function  of  time  in  Figure  9.  The  following  beam 
parameters  were  used:  beam  injection  energy  = 511  keV,  charge  per  simu- 
lation particle  = 70.78  statcoulombs , radial  position  for  particle  in- 
jection = 2.29  and  2.78  cm,  Z spacing  per  particle  = 0.25  cm,  and  beam 
current  risetime  0.15  ns.  At  the  beginning  and  end  of  the  pulse  the  Z 
spacing  of  particles  is  stretched  out  to  provide  an  approximate  exponen- 
tial rise  and  fall  of  the  beam  current.  The  saturation  beam  current  is 
4900  amps  that  provides  an  impedance  match  to  the  inner  transmission 
line  at  the  accelerating  voltage  of  300  kV. 

Two  problems  arise  when  using  the  code  with  a beam  of  charged  par- 
ticles. The  first  problem  appears  where  the  beam  is  injected.  Since 
the  code  formulation  assumes  continuity  of  charge,  an  injection  of  a 
charged  beam  implies  that  opposite  charge  is  building  up  at  the  surface 
of  the  injection.  Although  the  fields  due  to  this  surface  charge  decay 
rapidly  spatially  as  the  distance  from  the  injection  surface  increases, 
the  fields  are  sufficently  strong  to  produce  a sizable  perturbation  of 
the  beam  energy.  The  possibility  of  injecting  the  beam  "clothed"  with 
fields  it  would  have  traveling  in  an  infinite  pipe  was  considered,  but  a 
more  easily  implemented  solution  was  finally  used.  The  solution  is  to 
keep  the  particles  inactive  for  more  than  just  the  first  half  grid  spac- 
ing; they  are  kept  inactive  for  15  cm  in  our  case.  As  mentioned  pre- 
viously an  inactive  particle  is  one  that  produces  electromagnetic  fields 
but  does  not  respond  to  them.  This  procedure  allows  the  beam  to  pass 
through  the  region  of  unwanted  fields  undisturbed.  The  only  drawback  is 
that  the  diagnostic  that  checks  energy  conservation  must  be  modified  to 
account  for  the  "input  energy"  of  non-el ectromagnetic  origin  responsible 
for  maintaining  constant  velocity  of  the  beam  when  passing  through  the 
injection  surface  fields.  Such  a modification  has  not  been  implemented 
yet. 


Vacuum  Tube  Amplifiers,  George  E.  Valley , Jr. , and  Henry  Wallman , p. 
66  and  p.  77-78 , McGraw-Hill,  (1948). 
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Figure  8.  BQ  6 ns  after  switch  closing  in  the  baffled  cavity.  Every  other  point  is 
plotted  in  the  Z direction. 


2.7- 


Fiqure  9.  Energy  of  the  accelerated  beam  in  units  of  me  as  a function 
of  time  after  switch  closing.  The  injected  beam  has  . = 2. 
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The  second  problem  related  to  simulation  with  a beam  is  numeric: 

13  14 

the  Cherenkov  instability.  ’ Basically  this  instability  is  due  to 
the  dispersion  properties  of  the  finite  difference  equations.  Since  the 
higher  frequencies  propagate  at  less  than  the  speed  of  light,  it  is  pos- 
sible for  particles  traveling  at  less  than  the  speed  of  light  to  travel 
faster  than  those  higher  frequencies.  When  this  happens  the  particles 
radiate  into  these  higher  order  modes  via  the  Cherenkov  process.  This 
radiation  can  reach  very  high  levels  and  can  result  in  huge  perturbations 
to  the  particle  energies.  In  our  case  this  radiation  rises  rapidly  from 
, about  60  cm  from  the  beam  injection  surface.  We  have  been  able  to  live 

with  this  problem  by  diagnosing  the  accelerated  beam  energy  at  40  cm 
from  the  injection  surface  at  which  point  the  Cherenkov  instability  is 
of  negligible  magnitude.  Nevertheless  this  instability  is  a serious  ob- 
stacle to  the  investigation  of  a beam  traveling  through  several  cavities 
in  series.  It  seems  that  an  entirely  different  formulation  might  be  re- 
quired for  such  a case. 


V.  DISCUSSION 

The  above  results  have  shown  that  the  gross  features  of  the  simula- 
tion code  agree  rather  well  with  the  idealized  transmission  line  analysis 
for  the  case  of  the  3 to  1 impedance  ratio  investigated.  In  particular 
the  simulation  code  shows  how  well  the  pulses  turn  the  various  corners 
involved.  The  simulations  have  shown  that  the  generation  of  physical 
overshoots  and  higher  order  mode  excitation  depends  rather  sensitively  on 
the  detailed  geometry  of  the  switch.  For  example,  movement  of  the  lower 
boundary  of  the  switch  one  radial  position  resulted  in  significantly 
worse  overshoots  of  the  open  circuit  voltage. 

Calculations  have  been  performed  for  the  effects  of  the  discontinu- 
ities on  the  risetimes  of  the  pulses  by  employing  equivalent  circuit 

techniques.  Such  an  analysis  leads  to  a rise  that  is  not  a single  ex- 
ponential but  is  a sum  of  exponentials.  Moreover  at  the  beginning  of 
the  pulse,  a sharp  spike  in  the  reverse  direction  appears.  In  our  simu- 
lations no  such  reverse  spike  occurs,  although  the  simulation  is  consis- 
tent with  the  other  features  of  the  equivalent  circuit  analysis.  One  is 
therefore  led  to  the  definite  conclusion  that  the  reverse  spikes  are  ar- 
tifacts of  the  equivalent  circuit  model  and  the  model  does  not  properly 
predict  very  high  frequency  behavior  of  the  physical  system. 

The  importance  of  synchronization  of  the  beam  pulse  with  the  switch 
closing  was  mentioned  earlier.  If  the  beam  is  injected  too  early  the 
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dl/di  induced  voltage  will  decelerate  the  leading  portion  before  the  ac- 
celerating pulse  has  arrived  at  the  gap.  And  if  the  beam  arrives  too 
late,  then  the  leading  portion  of  the  beam  will  see  the  full  open  cir- 
cuit accelerating  voltage  and  receive  too  much  acceleration.  A differ- 
ence of  only  0.3  ns  in  the  timing  of  the  beam  injection  is  sufficient  to 
produce  significantly  poorer  results  than  shown  in  Figure  9.  For  many 
applications  a constant  accelerated  beam  energy  is  required.  If  magnetic 
steering  is  employed  in  the  beam  transport  system  for  example,  then  such 
a simple  consideration  as  reduction  of  extraneous  radiation  dictates  the 
need  for  a constant  energy  beam.  Thus  anything  affecting  the  constant 
beam  output  energy  is  a critical  design  issue. 

At  this  point  it  is  in  order  to  discuss  the  various  limitations  of 
our  computer  simulations,  their  motivation  and  justifications.  The  most 
apparent  limitation  is  the  two-dimensional  nature  of  the  code.  The  moti- 
vation for  this  limitation  is  simple;  a 3D  code  would  require  excessive 
computer  resources  at  the  present  time.  Also  many  interesting  geometries 
can  be  approximated  as  2D  structures  and  be  analyzed  economically  with  a 
2D  code.  The  geometry  analyzed  in  the  present  report  could  in  practice 
be  made  rotationally  symmetric  except  for  the  switches;  and  if  several 
switches  were  used,  they  should  closely  approximate  a 2D  switch  such  as 
that  assumed  here. 

Another  limitation  is  the  consideration  of  only  TM  modes.  Since 
there  is  nothing  in  the  geometry  to  excite  transverse  electric  (TE) 
modes  or  to  couple  TE  and  TM  modes,  it  would  appear  wasteful  of  computer 
resources  to  include  the  TE  modes  as  including  them  would  complicate  the 
programming,  nearly  double  the  storage  requirements,  and  more  than  double 
the  running  times  of  the  program.  It  must  be  stated, however, that  in  a 
working  accelerator  a stabilizing  axial  magnetic  field  would  be  neces- 
sary and  would  impart  helical  motion  to  the  particle  beam.  The  leading 
and  trailing  ends  of  a beam  with  this  motion  would  excite  some  TE  modes. 
Such  excitation  is  expected  to  be  very  weak  and  not  of  significance. 

Together  with  all  codes  using  finite  difference  techniques  for 
solving  partial  differential  equations,  the  present  code  requires  a fin- 
ite grid  spacing  and  finite  time  step  and  these  limit  the  spatial  and 
temporal  resolution  of  any  simulation.  Also,  due  to  the  resulting  dis- 
persion properties  of  the  finite  difference  equation,  the  accuracy  of 
the  simulation  is  less  than  perfect;  and  for  the  fine  detail  errors  of 
10%  can  be  expected.  The  finite  temporal  resolution  prevents  the  simula- 
tion of  very  high  frequency  instabilities,  and  other  techniques  are  re- 
quired to  verify  that  such  instabilities  can  be  avoided  in  a working  de- 
sign. Fortunately  analytical  procedures  exist  for  investigating  some 
such  possible  instabilities. 

The  final  limitations  of  the  simulations  that  we  discuss  are  those 
related  to  the  particle  beam.  The  main  simplification  made  to  the  par- 
ticle motion  is  to  allow  only  axial  motion.  Any  helical  motion  is  neg- 
lected because, as  mentioned  above, it  is  expected  to  be  small  and  not 
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have  significant  effects.  Not  allowing  radial  motion  of  the  particles  is 
motivated  by  the  desire  not  to  simulate  any  axial  magnetic  field  at  this 
point.  That  is  to  say,  if  the  particles  had  a radial  degree  of  freedom, 
an  axial  magnetic  field  would  be  required  to  prevent  radial  space 
charge  blow-up  of  the  beam.  Additional  numeric  problems  may  be  intro- 
duced under  these  conditions  and  it  was  decided  to  perform  the  analysis 
without  these  complications  at  this  point  and  to  defer  investigation  phe- 
nomena due  to  radial  particle  motion  to  a later  study.  Another  limita- 
tion is  that  the  numeric  Cherenkov  instability,  as  discussed  earlier, 
poses  a severe  restriction  to  the  analysis  of  high  current  beams.  It  has 
been  determined  empirically  that  for  low  current  beams  the  build-up  time 
for  this  instability  increases  considerably, thereby  allowing  useful  simu- 
lation of  such  beams  over  much  longer  distances.  It  is  expected  that  the 
numeric  Cherenkov  instability  would  become  worse  if  radial  motion  of  the 
beam  were  permitted. 


VI.  CONCLUSIONS 

Even  though  the  present  simulation  code  has  the  various  limitations 
just  mentioned,  the  results  of  the  present  simulations  are  interesting 
and  useful.  They  have  shown  the  sensitivity  of  the  results  to  the 
switch  geometry,  provided  quantitative  results  on  the  risetimes  to  be 
expected  from  this  system,  and  shown  that  the  beam  couples  resistively 
to  the  accelerating  structure.  The  direction  for  future  work  is  clear: 
generalization  of  the  code  to  remove  its  limitations.  Generalization  of 
the  boundaries  to  allow  for  boundaries  not  along  the  coordinate  axis 
and  the  allowing  of  radial  motion  to  the  beam  should  be  the  first  items 
undertaken,  after  which  new  geometries  such  as  radial  pulse  lines  should 
be  investigated.  Various  approaches  to  the  numeric  Cherenkov  instability 
problem  should  be  investigated  to  determine  the  most  suitable  one  to  ap- 
ply to  our  cases.  Finally  certain  interesting  geometries  are  3D  in  na- 
ture and  such  problems  should  be  investigated  if  the  necessary  resources 
become  available. 
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APPENDIX 


PROGRAM  USERS  GUIDE 


In  this  appendix  we  give  a description  of  the  overall  code  logic, 
the  input  variables,  and  the  major  variables  used  in  the  code.  In  addi- 
tion the  control  cards  required  to  run  the  code  on  the  BRL  Cyber  76  will 
be  shown.  Although  a listing  of  the  Cyber  Control  Language  (CCL)  pro- 
cedure and  the  in-house  coded  Poisson  solver  will  be  represented,  the 
entire  code  listing  will  not  be  included  since  the  code  was  not  released 
to  us  from  NRL  with  the  intent  of  general  distribution. 

A.  Subroutine  Structure 


Secondary 

Subroutine  Entry  Points 


Function 


PROP 


Main  calling  routine  to  allow  for  vari- 
able dimensions  with  the  use  of  CCL 


i 


PROJEC 


RELAX 

PROPIN 


PUSHIN 

PPRINT 
CP  LOT 

SORCIN 


PROBEB 

DIAGEB 

REDSEB 

PUSHP 


SOURCE 


Actual  main  routine,  handling  most  of 
the  input,  setting  up  the  masks  and 
calling  other  routines 

Poisson  solver 

Initializes  electromagnetic  field  in- 
tegrator 

Integrates  the  electromagnetic  fields 

Diagnoses  electromagnetic  fields  by 
calculating  energies  and  calling  list- 
ing and  plotting  routines 

Not  used 

Initializes  particle  pusher 
Integrates  particle  motion 
Field  listing 

2D  field  intensity  plotter  for  the 
line  printer 

Initializes  particle  injection  / 

Injects  particles  j»‘ 
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Subroutine 


I unct i on 


Second;!  ry 
I int ry  I’o  i nts 


niAi’.PR 

Not  used 

ZERO 

Zeros  fields 

INJECT 

Not  used 

B.  Input  Data 

All  of  the  input  to  the  program  (except  some  dimension  information 
via  CCL)  is  with  NAMELIST  read  statements.  The  program  contains  four 
namelists:  NLMAIN,  NLPART,  NLMASK,  and  NLEBPR.  The  first  three  are  in 
PROJEC  and  the  last  one  is  in  PROPIN.  Each  namelist  is  read  once  except 
for  NLMASK  that  is  read  twice,  before  and  after  the  call  to  RELAX. 


NAMELIST  NLMAIN: 
Initialized 


Variable 

Value 

Function 

NF 

7 

Number  of  fields  plus  masks 

DR 

1 

R grid  spacing 

DZ 

1 

Z grid  spacing 

DT 

.25 

Time  steps 

C 

1 

Velocity  of  light 

K0UT 

5 

Number  of  major  loops  - field  listing 
and  plottings  may  be  given  at  the  end 
of  each  major  loop 

KINN 

1 

Number  of  minor  loops  - the  number  of 
time  steps  per  major  loop 

KFDI 

1 

Diagnostic  printer  control 

KFPR 

0 

Not  used 

KFRE 

0 

Not  used 

KPDI 

0 

Not  used 

KPPR 

0 

Not  used 

KPRE 

0 

Not  used 

Variable 


Initial i zed 
Value 


Function 


POT (I) 

0 

An  array  that  sets  the  initial  electro- 
static potential  on  the  conductors. 

The  subscript  corresponds  to  the  sec- 
ond subscript  of  NBOIJN1)  (see  below  in 

NLMASK) 

LPOT 

F 

If  true  call  RELAX 

LPC 

F 

If  true  plot  initial  potential 

LPP 

F 

If  true  list  initial  potential 

LERC 

F 

If  true  plot  initial  ER  field 

LERP 

F 

If  true  list  initial  ER  field  \ 

LEZC 

F 

If  true  plot  initial  EZ  field 

LEZP 

F 

If  true  list  initial  EZ  field  • 

LINT 

F 

If  true  integrate  field  along  Z direc- 
tion to  obtain  accelerating  potential 

NIL 

1 

Innermost  radial  index  for  LINT  inte- 
gration 

NIH 

1 

Outermost  radial  index  for  LINT  inte- 
gration 

NIP 

1 

Particle  energy  diagnosed  or  LINT  in- 
tegration every  NIP  time  steps 

NAMELIST  NLP ART : 

NPMAX 

During  proce- 
dure call 

Don't  modify  with  a read! 

NQMAX 

During  proce- 
dure call 

Don't  modify  with  a read! 

NPS 

3 

Number  of  storage  locations  requried 
per  particle 

NZCEL 

8 

Number  of  particles  injected  per  Z , 

grid  spacing 

NRCEL 

8 

«• 

Number  of  particles  injected  per  R 
grid  spacing 
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Variable 


Initial i zed 
Value 


Function 


Ml  NR 

10 

Minimum  radius  of  injected  particles 

MAXR 

15 

Maximum  radius  of  injected  particles 

NRl’S 

0 

Number  of  different  particle  radii  -- 
redefined  after  the  read 

RPOS  ( 1 ) 

0 

Radii  of  the  particles  - redefined 
after  the  read 

GAMMA 

2 

Injected  particle  y 

SVELO 

1 

Sign  of  injection  velocity 

Wl’LAS 

1 

Plasma  frequency  - redefined  if  PCHAR 
^0 

PMASS 

1 

Particle  mass 

SCHAR 

-1 

Sign  of  particle  charge 

QCHAR 

Particle  charge  - output  only 

PCHAR 

0 

Particle  charge  - if  not  zero  used  as 
particle  charge  otherwise  calculated 
from  WPLAS 

PPAR1 

0 

Rise  and  fall  time  of  the  injected 
beam  in  time  steps 

PPAR2 

0 

Not  used 

PPAR3 

0 

Not  used 

PPAR4 

0 

Not  used 

PPAR5 

0 

Not  used 

PPAR6 

0 

Not  used 

PPAR7 

0 

Not  used 

PPAR8 

0 

Not  used 

JSTART 

0 

Time  step  to  start  injecting  particles 

JSTOP 

0 

Time  step  to  stop  injecting  particles 
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Variable 


Initial! zed 
Value 


Function 


LPV 

F 

If  true  diagnoses  particle  energies 

IPVU 

1 

Z position  to  diagnose  particle  ener- 
gies 

NAMELIST  NLMASK : 

NZMASK 

F 

If  true  mask  EZ  field 

NRMASK 

F 

If  true  mask  ER  field 

NTMASK 

Not  used 

NAREAS 

0 

Number  of  square  masks 

N BOUND ( I , J) 

0 

Corners  of  mask  in  grid  units  minus  1. 

Order  of  I-Z  . , Z , R . , R 

min  max  mm  max 

NAMELIST  NLEBPR: 

LBNZ 

F 

If  true  use  periodic  boundary  condi- 
tions 

IPRT 

6 

Unit  number  for  writing  energy  diagnos- 
tics 

1PLT 

6 

Not  used 

ISAV 

6 

Unit  for  storage  of  energy  diagnostics 

I FPR 

Not  used 

ISTR 

6 

Not  used 

I STW 

6 

Not  used 

NCFR 

1 

Frequency  of  plotting  in  R direction 

NCFZ 

1 

Frequency  of  plotting  in  Z direction 

NPFR 

1 

Frequency  of  listing  in  R direction 

NPFZ 

1 

Frequency  of  listing  in  Z direction 

NPWA 

1 

Not  used 

LSFR 

F 

Not  used 
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Initial!  zed 


Variable 

Value 

Funct i on 

LSPL 

F 

Not  used 

LPRT 

F 

If  true  print  energy  diagnostics 

I.SAV 

F 

Not  used 

LPLT 

F 

If  true  print  out  fields 

LFPR 

If  true  plot  fields 

LSTR 

F 

Not  used 

LF-- 

F 

If  true  plot  field  component  -- 

LP-- 

F 

If  true  print  field  component  -- 

NSW 

0 

Number  of  slow  risetime  switches 

RTIME 

1 

Risetime  of  the  switches  in  time  stej 
units 

ISW(I.J) 

0 

Switch  boundaries  (as  in  NBOIJND) 

C.  Major 

Variables  Not  Appearing 

as  I nput 

In  subroutine  PRO.JEC  all  of  the  fields  are  stored  in  the  array 
FL(I,J,K)  where  the  first  subscript  corresponds  to  the  Z position,  the 
second  to  the  R position,  and  the  third  to  the  field  type.  In  the  other 
subroutines  the  fields  are  stored  under  two  dimensional  arrays  (but  in 
the  same  core  locations)  EZ,  ER,  BT,  MEZ,  MER,  .JZ  and  RO;  where  BT  is 
B , MEZ  and  MER  arc  the  E and  E masks,  RO  is  p,  and  the  other  arrays 

have  obvious  names.  In  subroutine  RELAX  the  potential,  P,  is  stored  in 
the  core  locations  used  elsewhere  for  B^  thus  the  potential  is  not  avail- 
able after  the  electromagnetic  propagation  calculation  begins. 

The  array  CF(I,.J)  stores  the  radial  differencing  coefficients.  In 
RELAX  these  are  the  coefficients  in  eq.  (38);  and  for  the  propagation, 
the  coefficients  are  those  in  eq.  (19).  The  array  RPOS  contains  the 
radial  positions  of  the  particles  spaced  to  provide  a uniform  current 
density.  Information  of  the  particles  is  stored  in  the  array  PT(l)  for 
the  active  particles  and  in  the  array  QT(I)  for  the  inactive  particles. 

In  subroutine  PUSHIN  these  arrays  are  called  PSTR  and  QSTR . Each  three 
consecutive  elements  of  these  arrays  contain  the  Z position,  R position, 
and  specific  momentum  (momentum  per  unit  mass)  of  a particle.  As  des- 
cribed in  the  text,  the  particle  momentum  is  entirely  in  the  Z direction. 
The  variable  TIME  is  an  integer  variable  that  counts  the  number  of  time 
steps . 
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I).  Logic  Flow 


We  now  give  a brief  overview  of  the  logic  flow  of  the  program.  The 
main  routine  is  PROP  with  the  dimensions  for  the  fields  and  particle 
storage  assigned  from  the  CCL  procedure  call.  The  actual  main  driver, 
subroutine  PROJEC  is  called  and  first  initializes  the  NAMELIST  variable. 
Then  the  plasma  frequencies  or  particle  charges  arc  calculated  as  well 
as  the  RPOS  array.  The  mask  for  the  potential  solver  are  set  up  and  the 
potential  solver  is  called  if  requested.  Then  the  propagation  masks  are 
set  up;  and  the  initializing  portions  of  the  particle  pushing  E M pro- 
pagating, and  particle  subroutines  are  called.  Within  the  main  program 
loop  that  is  executed  next,  are  calls  to  the  F.+M  propagation  routine, 
the  particle  injector,  and  the  particle  pusher  as  well  as  diagnostics 
for  the  accelerating  potential,  beam  energy  and  field  values.  Liberal 
use  of  comment  cards  makes  the  program  flow  easy  to  follow. 

E.  Usage 

The  program  is  accessed  by  attaching  an  library  film  containing  a 
CCL  procedures  and  the  program  object  routines.  The  procedure  is  called 
by  its  name  followed  by  several  parameters  as  follows: 


Parameter 

1st/ 2nd 
Default 

Function 

NZ 

180 

Number  of  mesh  line  in  the  Z direction 

NR 

23 

Number  of  mesh  lines  in  the  R direction 

NPMAX 

45000 

Maximum  storage  for  active  particles 

NQMAX 

10000 

Maximum  storage  for  inactive  particles 

LL 

5000 

Print  limit  (lines) 

COMP 

F/T 

If  specified  recompile  one  or  more  subrou- 
tines 

LIST1 

L=0/SL, R=3 

List  the  main  routine  if  specified 

LIST2 

L=0/SL,R=3 

List  any  recompiled  routines  if  specified 

If  any  parameter  is  left  blank  in  the  procedure  call,  its  first  default 
is  chosen;  if  its  name  is  mentioned,  the  second  default  is  chosen;  and 
if  a numeric  value  or  a literal  inclosed  in  $ signs  is  used,  these  values 
are  substituted  for  the  parameter  in  the  body  of  the  procedure.  When 
COMP  is  specified,  the  procedure  itself  attaches  a file  containing  an 
UPDATE  library  of  the  source  routines  and  calls  UPDATE.  Only  those  rou- 
tines modified  are  recompiled;  the  others  are  loaded  from  the  object 
library.  When  recompiling,  the  section  following  the  control  cards 
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contains  the  UI’DATH  directives  followed  by  the  section  of  data.  Other- 
wise the  data  section  follows  the  control  card  section.  When  the  acccl 
crating  voltage  or  the  final  beam  y are  computed  and  written  to  TAI’liS 
or  TAPL9  respectively,  control  cards  should  be  included  for  the  disposi 
tion  of  these  files.  For  example,  they  may  be  cataloged,  copied  to  out 
put,  and/or  handled  by  calcomp  routines.  A listing  of  the  driving  pro- 
cedure, subroutine  RliLAX,  and  sample  inputs  follow.  The  sample  input 
shown  is  for  a special  diagnostic  run  with  the  baffle  geometry  shown  in 
Figure  7 and  does  not  initiate  beam  injection.  It  is  presented  merely 
to  show  the  overall  structure  of  the  input. 
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LIST  l=Y.L  = 0*/lSL.R=3l,.LlblZ  = S>L  = 0,i>/i>i>L.R=J4.WbMHbH  = \DATA. 

*■  TN » I =bribHHb .OPT  .LIST  1 .bYSEDIT.PL=LL. 

1 F t .COMP • CmnG. 

a t TACH.OLOPL.NHLSOURCt. 1U=ARL.PW  = *****. 

UPDATE  . 

f TN.I fOPT .LIST2.H=OUTtbYSEUlT. 

&E TURN. ULOPL .COMPILE. 

Endif.chng, 

MAP. PART. 

1FE.COMP.C2. 

luset.preset=nginf. 

LOAD. OUT . 

LOO. 

ELSF.C?. 

LOSE! .PRESt  T = NOINF. 
lGO. 

ENDIF ,C2. 

.DATA 

PROGRAM  PROP  ( INPUT .OUTPUT. TAPEH.T APER.TAPE5=INPUI . 

. TAPE<5  = OUTPUr  .UtbuO=UUTPUT) 
olMl-WSlON  FL  (NZ  *Nk  . 7 ) .CF  (NR. 2) 

dimension  pt (npmaxi .ot <nqmax> 

level  2.FL 
LeveL2.  PT.OT 
CJ  iMUN  PT.QT 
CO«lMUN  /FLL/  FU 

CAlL  PRO JtC  (FLtCF.NZ.NR. PT.OT .NPMAX, NoMAX) 

STOP  777 
END 
t 
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r.  c.  r. 


IBIS  PACE  IS  BEST  QUALITY  PRACTICABLE 
UK)!  OQPY  PURJUSHED  TO  DDC  


SUbKOUTINE  RELAX  INZ  .I.R.UZ  .OR.  ROT  *ER  » t Z . H » MSK.NA.Nb. 
. LRi-.LPP.LERC.LERR.  LtZC.LEZP 
. .KO.CA.Cb) 

DIMtNblON  EZ(NZ.NK).  tR(NZ.NR).  P(NZ.NN).  MSK(NZ.NR) 

. .ROKZO).  Nb(4.20) 

. .KO(NZ.NR).  CA(NK).  CH(nR) 

LEVEL  Z.  tW.tZ.P.MSN 

• . KO 

LOulLAL  L APR  lb)  . L Ah/  (b>  .1  Abl  (10)  .LAhC<9)  «LOC  (9)  .Lb  I (2) 
LOolLAL  lRC.LPR.  LHb,  LFKC.LERR.  LtZC.LEZP 
EOlil  VALE-ImCE  (LBC  ( 1 ) *L  ABC  ( 1 ) ) 

UAl.i  LI  T/HCPLU'>,"RK  IN"/ 

UA  I A L/ibR/”  K 1 •■ » Mi)I  Rt  . "C  T 1 0"  . “N  »•'/ 

l ) A 1 ii  LApZ/"  Z ."UIRK  "."Cl  10"  ."ix  -X."*"  I"/ 

0 A I A LAbC/"  "."l  li I , ••  AONO" ."ST  1C"."  E'OH". 

. « THt»,"  ” . ”f  I EL".  "0  !'•/ 

JAlrt  Lah  I EK  tT  ••.*•  tZ  HR  bl  BZ  •*. 

. 11  JR  •••••  JT  »,'•  JZ  «.'•  KO  ••/ 

data  1 1 R.Lbb  /ioh;  ; ( ( ; ; i itt » zh  p/ 

NR=0 


OD  2 0 j=NRMIN»NKMAX 
l.O  20  1 s.vZMIN.NZMAX 
NR=  IP* 1 

P( 1 .J)  = POT (N) 

^n  ms kii.j)  = o 

PM  = 0. 

00  Zb  1=1. NA 


RELAX 
KtLAX 
RELAX 
RtL  AX 
RELAX 
RELAX 
RELAX 

relax 

RELAX 

relax 

RELAX 

RELAX 

RELAX 

RELAX 

RELAX 

RELAX 

RELAX 

RELAX 

RELAX 

RELAX 


uo  10  J=1*NR 

RELAX 

Z2 

mi  10  1=1. NZ 

RELAX 

23 

tR ( i . J)  = 0. 

RELAX 

2a 

E / (1 . J ) = 0 , 

RELAX 

Zb 

R ( 1 . J)  = 0. 

relax 

Zb 

00  Z ' N=1 .NA 

relax 

27 

X/M  1 .<=  Nrl  ( 1 »N) 

♦ 1 

RELAX 

Zb 

N/«14X  = N t <2«N) 

♦ 1 

RELAX 

29 

NR'>'1N=  Nltl3.Nl 

♦ 1 

RELAX 

30 

NRMAX=  Nb(4.N) 

♦ 1 

RELAX 

31 

IF  (N/K1N  .LT. 

1)  N/min=1 

RELAX 

32 

IF  lIxZNAX  .01. 

NZ)  N/R.AX=NZ 

RELAX 

33 

IK  ( Nr Ml N .LT. 

1)  NrM 1 N= 1 

RELAX 

34 

IK  (NRMAX  .OT . 

NR)  NKKAX=NR 

RELAX 

3b 

RELAX 

RELAX 

RELAX 

RELAX 

RELAX 

RELAX 

RELAX 


25  PM 

= AMAXl  { PM .Abb  ( ROT ( I ) ) ) 

RELAX 

43 

IF 

(PM  .tU.  0.)  RM= 1 • 

RELAX 

44 

PI 

= Z.°  ASINIl.) 

RtL  AX 

4b 

Cl 

= 1.  / (2./0Z**Z  ♦ 2./UR«*2) 

RELAX 

46 

CC 

= L l/oZ**2 

RELAX 

4 / 

CO 

= (,:l/UR»«2 

RELAX 

4b 

Cc 

=—..ORI«Cl 

relax 

49 

C? 

= 1./  <2./UZ»"2  ♦4./DR»*2) 

RELAX 

bO 

CF 

= L2/b>Z**Z 

relax 

51 

CO 

= 4 • «C2/UR**2 

RELAX 

b2 

CM 

=-4.*RI»C2 

RELAX 

53 

00 

cb  1 =2.NR 

RELAX 

64 

CA ( 1 ) = (I-  .5)  *COZ  (1-1) 

RELAX 

55 

2b  CBU)  = (1-  1.5)  *CI)/  ( I - 1 ) 

RELAX 

bb 

RELAX 

57 

MAIN 

I Tt  RATION  SOLUTION  OF  R 

RELAX 

bB 

RELAX 

59 

OM 

= 1.7 

RtL  AX 

bO 

N7M1  =NZ-1 

RELAX 

61 

NRM 1 - NK-1 

RELAX 

62 

NP 

= INK-2)  * (N/-2)  - NP 

RELAX 

63 

46 


nrr,  n r.  o 


THIS  PACE  IS  BUST  QUALITY  PRACTICABLE 
FROM  COPY  FURBISHED  TO  DDC  y^- 


DO  30  LUUP31’400 

KELAX 

64 

DELV30. 

HELAX 

bb 

DELAV=0. 

HELAX 

66 

1)0  <*0  132*NZM1 

HELAX 

67 

DO  40  J=1,NRM1 

HELAX 

bB 

M = *n0(  mSK(I,J) >1) 

HbLAX 

69 

IF  (M  .bo.  0)  GO  TO  40 

HELAX 

70 

IF  U . b Q • 1 ) GO  TO  50 

HELAX 

71 

ObL  3 CC*  (P(1»1.J)  ♦ P(l-l.J))  ♦ CA(J)*P(1,J 

♦1)  *Cb(J)*P(I,j-|) 

HELAX 

72 

. ♦Cb»KU(l ,J)  - P(I*J) 

HELAX 

73 

GO  ro  60 

HELAX 

74 

So 

DEL  = CF*  (P(I*1.1)  ♦ P(l-l.l))  ♦ CG*P ( I • 2 ) 

♦ Ch*hO (1,1) 

HELAX 

75 

. -P< 1 .1) 

HELAX 

76 

60 

A'JtL  3 AH  SIDED 

HELAX 

77 

DbLV  = A MAX  1 (OELV  ,ADbL  > 

HELAX 

78 

Ob  LAV  = OELAV  ♦ AUEL 

HbLAX 

79 

PU.J)  = P(I*J)  ♦ DbL*OM 

HELAX 

bO 

40 

CONTINUE 

HELAX 

81 

Db'LV  = UELV/PM 

HELAX 

82 

obLAV  = DELAV/NP/PM 

HELAX 

83 

IF  (LOOP  . EQ.  1)  WHITE  (6,70) 

LOOP ,UbLV,UtLAV 

HtLAX 

84 

IF  (d.OOP/10)  *10  .bo.  LOOP)  whIIE  ( b , 70 ) 

LOOP, ObLV , OELAV 

HELAX 

85 

70 

FO«MAI  (••  LOOP  NO."  , 1 4 , ••  LAHGbST  bPb"* 

lPblb.4, 

hELAX 

86 

. ••  AVbHAGt  EPS",  b 1 5. 4 ) 

HELAX 

87 

IF  ( ObLV  .Lf.  l.b-7)  GO  TO  SO 

HELAX 

b8 

30 

CONTINUE 

HbLAX 

89 

HbLAX 

90 

CALCULATE  E FIELDS 

HELAX 

91 

HELAX 

92 

80 

DO  110  J=2.NK 

HELAX 

93 

bHU.J)  =<P(1.J-1)  - P(1,J>)  /OH 

HELAX 

94 

Du  110  I=2*NZ 

hELAX 

9b 

eru,j>  =(p(i«j-i)  - Pinjii  /or 

HELAX 

96 

no 

E7(i«J)  3<P(1-1,J)  -P(I»J))/DZ 

HELAX 

47 

DO  120  1=2»NZ 

RELAX 

98 

120 

EZU.l)  = (P(l-lill-  P(I,1))/DZ 

HELAX 

99 

HELAX 

100 

PAPER  PLOTS 

HELAX 

101 

HELAX 

102 

NE  = NZ/120* 1 

HELAX 

103 

ME  3 NH/100*1 

HELAX 

104 

LAHC ( 1 ) =LBT ( 1 ) 

HELAX 

105 

IF  (.NOT.  LPC)  GO  TO  90 

HELAX 

106 

LHC  I 7)  = LbB 

HELAX 

107 

CALL  CPLOT  ( P.NZ.lfNZ.NE,  NH , 1 » NH • ME » U.,0., 

LAbL,LAHZ,LAHH, 

HELAX 

108 

. M5K,  .b  .) 

HELAX 

109 

90 

IF  (.NOT.  LERC)  GO  TO  130 

HELAX 

no 

LrtC ( 7)  = LABI  (1) 

HbLAX 

111 

CALL  CPLOT  (EH,NZ  » 1 *NZ  « NE  » NH , 1 , NH • Mb  • U.,0., 

LAHCLABZ  ,labh. 

HELAX 

112 

. M6N , . 1 • ) 

HELAX 

ID 

130 

IF  (.not.  LEZC)  GO  TO  140 

HELAX 

n<* 

LHC ( / ) = LABI (3) 

hELAX 

lib 

CALL  CPLOT  <EZ*NZ*1 »NZ,NE»  NH , I , NH , ME , 0.,0.» 

L ABC ,LABZ,LAbH, 

HELAX 

lib 

. MSK,  .T.) 

RELAX 

117 

140 

COul 1NUE 

KELAX 

118 

LA-iC  ( 1 ) =LBT  (2) 

HELAX 

119 

lb  (.NOT.  LPP)  GO  TO  100 

RELAX 

120 

LHC ( n = LBB 

HEcAX 

121 

CALL  PPH1NT  ( P,NZ , 1 ,NZ , 1 , NH ♦ 1 , NH , 1 , OtLABC) 

RELAX 

122 

100 

IF  (.NOT.  LEHP)  GO  TO  150 

HELAX 

123 

LBCI7)  * LABI (1) 

RELAX 

124 

CALL  PPRINT  (ER,NZ,1,NZ,1,  NH.l.NH.i,  0 ,LABC) 

KELAX 

12b 

150 

IF ( .NOT . LEZP)  GO  TO  160 

RELAX 

126 

47 


r n r 


I 


/ 


LHCI7)  = LABK3) 

CALL  KKHINT  (tZ tNZ • 1 »NZ • 1 « NK.l.IMK.l.  0*LAbC> 

K1  Sf  f MASK 

160  DD  1 70  J=1 .N« 

!)<>  170  l = l»NZ 
171'  MSMI.J1  = 1TH 
Kt I UKN 
ENl) 


KELAX  U7 
KtLAX  12b 
KELAX  129 
KELAX  UO 
KELAX  1J1 
KELAX  U2 
KELAX  1 J 3 
KtLAX  1JA 
KtLAX  US 
KELAX  13b 


•1N.STMKZ.T100. 

ACCOUNT .*•***«.  SHNIDMAN.  H390  Xb8b9 
mTTACH.HL.NML  UbJtCT. Il)=AKL.Pta=*«******. 

LIHRARY .ML. 

VUPWl'P. • « 10. 10. lOOOO.CUMP. 

( 

*io  smak r 
«I  PPHlNT.AO 

00  BO  1=1. MM 
H = 1-1.5 
OU  HO  J=1.NN 
HO  A < J. I)  = A ( J» I ) *K 
7 

SNLNA1N 
KP0I=1 . 

IJT  = 1 .t-i  1 * 

C=3.E 1 0 . 

KlNN=t>UO«KOUT=l* 

LP01 =T  » 

POT = 1 000. . O.iO.iU.iU.iU.i  iooo..o..looo..o.,iooo.«o..luoo.»o.,iouo..o., 

1000. »0.t 

s 

iNt.PAR  I 
NZCEL=a.. 

NMCEL=2. 

PCHAk=-/O.770A, 

PMASS=l .3A2296E-lb, 

PPAW 1 = 1 b. . 

MAXW=3. 

M I Nk  = ? « 

JST  ART  = 10000. 

S 

tNLMAS* 

NZMASK=T . 

NRMASK= I . 

NARt As  = 1 h .NHOUND ( 1 . 1 ) = 2b . 1 7 1 . 1 a . 1 S . 2b«179*A*b»  170.1/9.0.22.  O.l.O.S* 

0.18.A.2. . lb. 179*21 .22. 

20.30.0. 1**.  33. 37. b. 11.  AU.AA.6.1A.  A7»bl*b»ll*  SA.Se.b.lA.  bl.bb.b.ll. 

68.72.0.  1A.  7b. 79. 5. 11.  b2.8b,b,lA.  b9«93.b.ll.  96.100.6.1a,  1 03» 107 .b. 1 1 . 

S 

SNLMAS* 

i > 

SNLEHPR 

LPWT=T. 

LPfc  R=T. 

LPLT=T . 

NCFZ=2. 

NSw= 1 • lbW=lb.27.lA,21. 

WTlME=lb.. 

S 

? 
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